* Gary LaMotte's Mathematica solutions to exercises from the 2012 Gravitational Wave Astronomy Summer School

*DAY 2*

*1.  PREPARATIONS*

*PARTS   a. and   b.*

samprate=16384;

template= Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day2template.txt", "List"];

templatesecs = Length[template]/samprate

ttemplate =  Table[(s-1.0)/samprate  , {s,Length[template]}];

templateseries = Transpose[{ttemplate, template}];

ListLinePlot[templateseries, PlotLabel -> Template, AxesLabel ->  {Time (Secs), Amplitude}, GridLines->Automatic]

*TEMPLATE HAS A STARTING FREQ OF APPOX 40Hz AND FINAL MERGER FREQ OF APPROX 500 Hz*

*PART  c.*

data =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day2data.txt", "Table"]];

ndata =Length[data]

datasecs = Length[data]/samprate

tdata =  Table[(r-1.0)/samprate  , {r,Length[data]}];

dataseries = Transpose[{tdata, data}];

ListLinePlot[dataseries, PlotLabel -> Data, AxesLabel ->  {Time (Secs), Amplitude}]

*NO SIGNAL SEEN*

*Part d.*

firstpartdataseries = dataseries[[1;;10000]];

ListLinePlot[firstpartdataseries, PlotLabel -> Data, AxesLabel ->  {Time (Secs), Amplitude}]

*Part e.*

*START WELCH METHOD ANALYSIS*

*Definitions to use to insert into Welch Method "Subprogram" are as follows:*

ts = data;  nsec = datasecs;

*Then plug into Subprogram*

partitions = 8;

segments = Partition[ts,     Length[ts]/partitions,      Length[ts]/(2*partitions)];

nsec2 = nsec / partitions

ndata2 =  (samprate * nsec) / partitions

newtime = Table[(p-1.0)*(nsec2/ndata2)  , {p,ndata2}];

filtsegs = Transpose[(0.54 - 0.46(Cos[2*Pi* newtime/nsec2]))*Transpose[segments]];

ftfiltsegs = (Abs[Fourier[filtsegs]])^2;

nyqtrans = Transpose[ftfiltsegs][[1;;(samprate/1)]];

nyqpart = Transpose[nyqtrans];

added= Total[nyqpart];

averaged = added/((partitions*2) - 1);

index = Table[(p-1.0)*partitions /datasecs, {p,Length[averaged]}];

normalized = Transpose[{index, averaged}];

ListLogPlot[normalized,PlotRange -> {{0,8192},{.001, 1000}}, Joined->True, PlotLabel-> PSD Welch Data, GridLines-> {{1000,2000,3000,4000,5000,6000,7000,8000}, {.01,.1,1,10,100,1000}}, AxesLabel->{Frequency, Power}]

ListLogPlot[normalized,PlotRange -> {{0,1024},{.01, 1000}}, Joined->True, PlotLabel-> PSD Welch Data]

*END WELCH*

*2.  A VERY SIMPLE DIGITAL FILTER*

*Part a.*

shift = Drop[data,-1];

shiftedbyone = Insert[shift,shift[[1]],1];

*BELOW IS THE SIMPLE  FIRST-DIFFERENCE  FILTER*

subfilt = data - shiftedbyone;

st = Transpose[{tdata, subfilt}];

*Part b.*

firstst = st[[1;;10000]];

ListLinePlot[firstst, PlotLabel -> Filter Difference, AxesLabel ->  {Time (Secs), Amplitude}]

*Very low frequency oscillation has been removed (10 to 20 Hz which was the dominant freq seen in the PSD in Part 1).*

*Part c.*

*Welch Method*

*Definitions to use to insert into Welch Method "Subprogram" are as follows:*

ts = subfilt;  nsec = datasecs;

partitions = 8;

segments = Partition[ts,     Length[ts]/partitions,      Length[ts]/(2*partitions)];

nsec2 = nsec / partitions

ndata2 =  (samprate * nsec) / partitions

newtime = Table[(p-1.0)*(nsec2/ndata2)  , {p,ndata2}];

filtsegs = Transpose[(0.54 - 0.46(Cos[2*Pi* newtime/nsec2]))*Transpose[segments]];

ftfiltsegs = (Abs[Fourier[filtsegs]])^2;

nyqtrans = Transpose[ftfiltsegs ][[1;;  Length[Transpose[ftfiltsegs ]]/2]];

nyqpart = Transpose[nyqtrans];

added= Total[nyqpart];

averaged = added/((partitions*2) - 1);

index = Table[(p-1.0)*partitions /datasecs, {p,Length[averaged]}];

normalized = Transpose[{index, averaged}];

ListLogPlot[normalized,PlotRange -> {{0,8192},{.0000001, 1000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch, GridLines-> {{1000,2000,3000,4000,5000,6000,7000,8000}, {.000001,.00001,.0001,.001,.01,.1,1,10,100,1000}}, AxesLabel->{Frequency, Power}]

*Diminshed low freq noise.  Slightly increased high freq noise versus unsubtracted data.*

ListLogPlot[normalized,PlotRange -> {{0,1024},{.0000001, 1000}}, Joined->True, PlotLabel-> SegmentedHamming PSD Welch]

*End Welch*

*3.  A BETTER DIGITAL FILTER*

*PART   a.*

*BELOW:  SIMULATED BUTTERWORTH FILTER BY USING SQUARE BANDPASS FILTER 35 TO 250Hz*

*SUCCESSFUL SQUARE WINDOW FILTER DESIGN IN FREQUENCY DOMAIN,INVERSE FOURIER TRANSFORMED TO THE TIME DOMAIN (SINC FUNCTION)AND CONVOLVED  WITH TIME SERIES*

*SQUARE FILTER PART*

*35Hz is mapped to 560 by multiplying by a factor of 16 and 250Hz is mapped to 4000 in the same manner.  This technique allows the inverse fourier transformed square bandpass filter to be convolved with the data.*

*The following are inputted into the Square Filter Subprogram:*

upcorn =   250;

lowcorn =  35;

*Square Filter Subprogram starts here.*

start = {1,1};

pass = PadRight[start, ((upcorn - lowcorn)*datasecs)+1, 1];

rtside = PadRight[pass, ndata-((lowcorn*datasecs)-1),0];

sqfilt = PadLeft[rtside,ndata];

ListLinePlot[sqfilt, PlotRange->{{0,5000},{0,2}}]

filtergraph = Transpose[{ samprate*tdata/datasecs,sqfilt}];

ListLinePlot[filtergraph, PlotRange->{{0,500},{0,2}}]

*END SQUARE SUBPROG.*

sinc = InverseFourier[sqfilt];

*Below:  1/512 is a normalization constant*

realsinc = (1./512.) * Re[sinc];

convts = ListConvolve[realsinc,data,1];

sincts = Transpose[{tdata, convts}];

*PART   b.*

ListLinePlot[sincts, PlotLabel -> Filter BandPass, AxesLabel ->  {Time (Secs), Amplitude},PlotRange->{{0,16}, {-.3,.3}}]

*A signal is seen with this square filter in Mma just like with the Butterworth filter using Matlab.*

*4.  MATCHED FILTERING IN THE TIME DOMAIN*

*PART   a.*

*First pad the template vector with zeros so it's equal in length to the data vector.*

templateoo = PadRight[template,262144];

*Then step the template through the data to obtain a correlation at each step.*

table = Table[templateoo . RotateLeft[data,i ],{i,0,(Length[data] - Length[template])}];

finaltable = PadRight[table,262144];

correlationseries = Transpose[{tdata, finaltable}];

*PART   b.*

ListPlot[correlationseries , Joined -> True, PlotRange-> {{0,16}, {-.00000000000000006,.00000000000000011}}, AxesOrigin->{-.00000000000000006,0},  Frame->True,FrameLabel->{Time,Magnitude}]

*Obvious strong signal starting just after 5 secs and a much smaller possible signal later just after 13 secs.   Note that the data in the last second is not evaluated as the template has halted it's stepped progression through the data.  The template is one second long.*

*5.  DOWN-SAMPLING*

*PART   a.*

StandardDeviation[data]

*PART   b.*

downsampfactor = 16;

sixteenth = Partition[data,1,downsampfactor];

downsamp = Flatten[sixteenth];

StandardDeviation[downsamp]

*PART   c.*

*The following are inputted into the Square Filter Subprogram:*

upcorn =   512;

lowcorn =  0;

*Square Filter Subprogram starts here.*

start = {1,1};

pass = PadRight[start, ((upcorn - lowcorn)*datasecs)+1, 1];

rtside = PadRight[pass, ndata-((lowcorn*datasecs)-1),0];

sqfilt = PadLeft[rtside,ndata];

ListLinePlot[sqfilt, PlotRange->{{0,10000},{0,2}}]

filtergraph = Transpose[{ samprate*tdata/datasecs,sqfilt}];

ListLinePlot[filtergraph, PlotRange->{{0,600},{0,2}}]

*END SQUARE SUBPROG*

sinc = InverseFourier[sqfilt];

*Below:  1/512 is a normalization constant*

realsinc = (1./512.) * Re[sinc];

convts = ListConvolve[realsinc,data,1];

sincts = Transpose[{tdata, convts}];

ListLinePlot[sincts, PlotLabel -> LowPass Filter, AxesLabel ->  {Time (Secs), Amplitude},  PlotRange->{{0,16},{-2,2}}]

sixteenthconvts = Partition[convts,1,downsampfactor];

downsamp2 = Flatten[sixteenthconvts];

tdownsamp2 =  Table[(r-1.0)*16/samprate  , {r,Length[downsamp2]}];

dataantialiasdownsamp = Transpose[{tdownsamp2, downsamp2}];

ListLinePlot[dataantialiasdownsamp, PlotLabel -> Filter Lowpass Downsampled, AxesLabel ->  {Time (Secs), Amplitude},  PlotRange->{{0,16},{-2,2}}]

*PART   d.*

StandardDeviation[downsamp2]

*So the Std Dev decreased after applying a low pass filter as expected.*

*PART   e.*

*Welch Method*

*Definitions to use to insert into Welch Method "Subprogram" are as follows:*

ts = downsamp2;  nsec = datasecs;

partitions = 8;

segments = Partition[ts,     Length[ts]/partitions,      Length[ts]/(2*partitions)];

nsec2 = (nsec/downsampfactor) / partitions;

ndata2 =  (samprate * (nsec/downsampfactor)) / partitions

newtime = Table[(p-1.0)*(nsec2/ndata2)  , {p,ndata2}];

filtsegs = Transpose[(0.54 - 0.46(Cos[2*Pi* newtime/nsec2]))*Transpose[segments]];

ftfiltsegs = (Abs[Fourier[filtsegs]])^2;

nyqtrans = Transpose[ftfiltsegs ][[1;;  Length[Transpose[ftfiltsegs ]]/2]];

nyqpart = Transpose[nyqtrans];

added= Total[nyqpart];

averaged = added/((partitions*2) - 1);

index = Table[(p-1.0)*partitions /datasecs, {p,Length[averaged]}];

normalized = Transpose[{index, averaged}];

ListLogPlot[normalized,PlotRange -> {{0,512},{.0000001, 1000}}, Joined->True, PlotLabel-> Downsamp Antialias PSD Welch, GridLines-> {{1000,2000,3000,4000,5000,6000,7000,8000}, {.000001,.00001,.0001,.001,.01,.1,1,10,100,1000}}, AxesLabel->{Frequency, Power}]

*End Welch*

*PART   f.*

*Downsample Template:*

downsamptemplate = Flatten[Partition[template,1,downsampfactor]];

tdownsamptemplate =  Table[(r-1.0)*16/samprate  , {r,Length[downsamptemplate]}];

*Matched Filtering*

downsamptemplateoo = PadRight[downsamptemplate,16384];

tabledownsamp = Table[downsamptemplateoo . RotateLeft[downsamp2,i ],{i,0,(Length[downsamp2] - Length[downsamptemplate])}];

finaltabledownsamp = PadRight[tabledownsamp,16384];

correlationseriesdownsamp = Transpose[{tdownsamp2, finaltabledownsamp}];

ListPlot[correlationseriesdownsamp, Joined -> True, PlotRange-> {{0,16}, {-.000000000000000002,.0000000000000000034}}, AxesOrigin->{-.00000000000000006,0},  Frame->True,FrameLabel->{Time,Magnitude}]

*Obvious strong signal starting just after 5 secs and a much smaller possible signal later just after 13 secs.   Note that the data in the last second is not evaluated as the template has halted it's stepped progression through the data.  The template is one second long.*